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Abstract: We present a method for high-resolution reconstruction of 
fluorescent images of the mouse thorax. It features an anatomically guided 
sampling method to retrospectively eliminate problematic data and a 
parallel Monte Carlo software package to compute the Jacobian matrix for 
the inverse problem. The proposed method was capable of resolving 
microliter- sized femtomole amount of quantum dot inclusions closely 
located in the middle of the mouse thorax. The reconstruction was verified 
against co-registered micro-CT data. Using the proposed method, the new 
system achieved significantly higher resolution and sensitivity compared to 
our previous system consisting of the same hardware. This method can be 
applied to any system utilizing similar imaging principles to improve 
imaging performance. 
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1. Introduction 

Fluorescence diffuse optical tomography (FDOT) is a noninvasive method for imaging 
fluorescent agents in living organisms using far-red or near infrared light. With the rapid 
development of fluorescently labeled molecules as the targeted imaging probes, FDOT 
promises great potential for preclinical studies [1-3]. Compared to other imaging modalities, 
such as magnetic resonance imaging (MRI), x-ray computed tomography (CT), and 
ultrasound, a major disadvantage of FDOT for preclinical imaging is the difficulty in 
modeling light propagation in highly scattering and heterogeneous media. Many techniques 
have been proposed to overcome this limitation by addressing one or more of the difficulties 
in modeling, e.g., the use of Born normalization [4], higher-order harmonics for the source 
term [5], mathematical approximation using simplified geometry [6], altered geometry of the 
animal using a compressing imaging chamber [7], and the time-gated early-photons [8]. 
Nonetheless, it is still a challenge for FDOT to resolve closely located small fluorescent 
heterogeneities in the torso of small animals, which is critical for effective preclinical 
imaging. The thorax of the small animal, in particular, is of great interest for preclinical 
applications because it contains the heart and lungs that are very important target organs for a 
large number of cardiovascular, pulmonary, and cancer disease models, which play critical 
roles in numerous preclinical and translational efforts. 

Light propagation in the turbid media, also known as the "forward problem," has been 
reported extensively in the literature [9-11]. The main category of the forward modeling 
methods is based on the radiative transfer equation (RTE). Due to the complexity of the RTE, 
the diffusion approximation is often applied in the context of diffuse optical imaging, 
including diffuse optical tomography (DOT) and FDOT. In practice, further simplification of 
the diffusion equation, most notably using the Born approximations, are frequently used to 
solve the forward problem in diffuse optical imaging [4]. Another important category of the 
forward models is based on the Monte Carlo method, which is known for its accuracy but also 
for its high computational cost [12-14]. Reconstruction of the FDOT image from a series of 
optical measurements and the knowledge of the optical medium, as well as the boundary 
conditions, is referred to as the "inverse problem." 

Successfully solving the inverse problem relies on solving the forward problem correctly, 
because the two problems share one component — light propagation in the medium. The 
specific choices of the forward and the inverse problem solvers define a particular 
reconstruction method for diffuse optical imaging. In preclinical imaging, the highly complex 
geometry typically requires the use of either the finite-element or Monte Carlo methods. On 
one hand, although the Monte Carlo method is computationally expensive compared to the 
finite-element method, it relies on fewer assumptions. For example, the diffusion 
approximation, on which most of the finite-element method is typically based, becomes 
invalid when the ratio of the absorption to the scattering coefficients is low (e.g., in the liver, 
blood, or bladder) or, when the source/detector distance is comparable to the transport mean 
free pathlength. In addition, the simplicity of the Monte Carlo method makes it attractive over 
the finite-element method, which involves complex boundary conditions. On the other hand, 
as the computer technology continuously improves, particularly the wide availability of 
parallel computing, the Monte Carlo method is becoming more and more attractive. In this 
study, the forward problem solver used a supercomputer-based parallel implementation of our 
Monte Carlo method [15]; and the inverse problem solver was based on the Born method. 

Our goal was to develop new data acquisition and image reconstruction techniques that 
can provide high spatial resolution and high sensitivity to satisfy the needs of preclinical 
applications. These techniques included (1) parallel implementation of our Monte Carlo 
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method that constructed the Jacobian matrix directly based on photon paths from the source to 
the detector, rather than taking the product of the two fields of photon density waves 
originated from the source and the detector separately; and (2) use of anatomy of the animal 
to automatically identify and eliminate problematic data points within a large data set. 

In diffuse optical imaging, the inverse problem is typically underdetermined in that there 
are usually fewer measurements than the number of unknowns. Therefore, for high-resolution 
reconstruction, denser spatial sampling is necessary. Using the free-space full-angle 
acquisition method [16], the inverse problem can become seemingly overdetermined using a 
scanning laser and a charge-coupled device (CCD) camera. Although studies show that higher 
sampling density can improve spatial resolution, e.g., in [17,18], sampling strategy plays an 
important role for efficient and high-resolution reconstruction [19-21]. Particularly in free- 
space animal imaging, where the exact surface of the animal is highly irregular and 
unpredictable, a mechanism is needed to accurately determine source/detector positions using 
an optical method during the experiment. We developed a geometrical calibration method to 
optically determine the source and detector positions and to subsequently co-register the 
optical and x-ray data. 

For diffuse optical imaging, the highly heterogeneous macroscopic anatomical structures 
of the animal significantly alter the propagation of light and thus heavily influence the 
detected signal. Consequently, incorporation of structural information is critical for correctly 
modeling light propagation and improving reconstruction in animal experiments, as concluded 
in many studies such as [22-26]. In this study, we not only used the x-ray micro-CT data 
registered to the optical data as the source of structural priors for FDOT reconstruction, but 
also used this information for optical data sampling. Specifically, instead of sampling 
uniformly, or using any predetermined pattern, we used an anatomically guided sampling 
strategy to exclude data that would otherwise inversely affect the reconstruction. 

2. Methods 

2.1. Animal preparation 

We chose quantum dots (QD) as the fluorescent contrast agent because of their high quantum 
yield, excellent photo stability, and broad selection of functional groups for targeted imaging 
in the future. Although the cytotoxicity of cadmium-based QD to living organisms is well- 
known, the use of long molecular chain polymer coatings, particularly polyethylene glycol, 
can effectively eliminate the toxicity concerns of the QD for in vivo studies, e.g., as reported 
in [27-29]. We used commercially available non-conjugated QD branded Qtracker 800 
(Invitrogen, Carlsbad, CA). At our excitation wavelength of 642 nm, the extinction coefficient 
of the QD is 750,000 M _1 cm _1 ; and the quantum yield is 3%, according to the product 
specifications. 

The animal experiment was conducted on a nude mouse (Nu/Nu, male, 25.5 g) cadaver to 
limit the complexity of animal support, while still providing a biologically realistic model. 
The QD inclusions were constructed by sealing polyethylene capillaries with paraffin and 
filled with 1 uL of 200 nmol/L QD saline solution, which were inserted into the thorax via the 
esophagus and the trachea from an incision on its neck. Each of these inclusions contained 
200 femtomoles of QD. The dimensions of the inclusions were 0.9 mm in a diameter and 1.6 
mm in length. The distance measured between their closest edges was 2.5 mm. The 
dimensions and the positions of the inclusions were verified by x-ray micro-CT. These 
inclusions were enclosed by the heart, the lungs, and the spine, representing a very challenge 
but nonetheless realistic situation for FDOT. The purpose of such a design was to emulate a 
reasonably favorable scenario for high-resolution FDOT study of the mouse with lung tumors: 
the dimensions of the inclusion are equivalent to a millimeter- sized tumor; and the same 
concentration of the inclusion is achievable in the blood of the mouse (a typical total volume 
of 2 mL) by injecting a reasonably large dose of QD (200 uL of 2 uM) via the tail vein. 
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2.2. Experimental setup 

Our optical imaging platform consists of a continuous-wave laser diode as the light source, a 
galvo-scanner that enables raster scanning of the laser on the animal, a rotation stage that 
vertically suspends the animal and provides unobstructed panoramic view of the animal, a 
CCD camera (Imager QE, LaVision, Germany) as a high-density photodetector array to 
measure the transillumination, and a filter wheel to select the appropriate detection 
wavelength for the camera. The control and synchronization of these devices were 
implemented using lab-developed Lab VIEW (National Instruments, Austin, TX) applications. 
This experimental setup is similar to the one previously described in [25], except that the 
diode laser and the filters were replaced to match the excitation and emission spectra of the 
QD. The laser diode (HL6388, Opnext, Fremont, CA) radiates at 642 nm with an optical 
output power of 250 mW. The absorption spectrum of the QD monotonically decreases from 
-300 nm in the ultraviolet band and extends beyond 700 nm in the near infrared band. The 
fluorescent emission spectrum of the QD peaks at -785 nm with a full- width at half- 
maximum bandwidth of -75 nm. Our choice of this particular excitation wavelength was 
based on considerations of optical penetration, the excitation efficiency, and the specifications 
of commercially available laser diodes. 

The experiment started with fluorescence acquisition, using two stacked 750 nm long-pass 
filters (CVI Melles Griot, Albuquerque, NM) to collect the fluorescent light and a 647 nm 
long-pass edge filter (Semrock, Rochester, NY) to block the laser. The CCD camera acquired 
one image, with integration time of 4 s, for every source position on the animal until the 
galvo-scanner completed a scan. The rotation stage subsequently changed the angular position 
of the animal, before the next laser scan and corresponding data acquisitions began. This was 
repeated until the rotation stage resumed its initial angular position. The same procedure was 
subsequently applied to acquire the excitation images using a neutral density filter to avoid 
CCD saturation (integration time 0.1 s). We used 6 vertically arranged source positions for 
each angular position, each separated by -2 mm and 40 angular positions, which resulted in 
240 images for both fluorescence and excitation. We applied time delays of 0.1 s between 
consecutive source points and 1 s between angular positions. Overall, fluorescence acquisition 
required 1024 s; acquisition of excitation images needed 88 s; and the geometrical calibration 
took 1.2 s. 

After the acquisition of fluorescence and excitation images, a set of reflectance images 
were taken for the same 40 angular positions using a diffusive LED white-light source 
(Thorlabs, Newton, NJ). These reflectance images were used to extract the 3-D surface of the 
animal using a parallel-beam backprojection algorithm. Details of this method were 
previously described in [25]. 

The last step in optical imaging was to acquire calibration images for precise spatial 
registration. A calibration target was constructed using a rigid translucent sheet material, on 
which a rectangular grid with known intervals was printed. The calibration target replaced the 
rotation stage and was positioned coplanar with the axis of rotation. The laser impinging on 
the target produced a bright spot. Under front illumination, both the printed grids and the laser 
spot were visible to the CCD camera. One CCD image was taken for each of the 6 
laser/gal vo- scanner positions. 

After optical acquisition, the rotation stage (with the animal undisturbed) was transferred 
to the adjacent x-ray system to acquire micro-CT data as structural information, which was 
then co-registered with the optical data using our calibration method. Details of our dual x-ray 
tube/detector micro-CT system were previously described in [30]. A set of 200 projections 
were acquired in less than 3 minutes using a single tube/detector with an angular step size of 
1.8° per projection. The x-ray tube operated at 80 kVp, 160 mA, and 10 ms per exposure. The 
radiation dose associated with the micro-CT scan was 7.2 cGy. The projection set was used to 
reconstruct the tomograms with a Feldkamp algorithm [31] implemented in the Cobra 
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EXXIM software package (EXXIM Computing, Livermore, CA). The micro-CT image was 
reconstructed as a (512) 3 matrix with an isotropic voxel size of 88 um. 

2.3. Data processing 

Based on the white-light reflectance CCD images, a sampling region of 9 x 18 mm was 
defined in the chest area of the animal in the emission and excitation CCD images. Within the 
sampling region, 4 x 10 sampling points were evenly spaced. The resulting sampling density 
was 3x2 mm horizontally and vertically, respectively. The mean value of the pixels within a 
lxl mm area centered at the sampling point was used as the optical signal intensity for the 
sampling point. The spatial sampling points were identical for the raw CCD images at the 
emission and the excitation wavelengths. The resulting excitation and emission data sets both 
had 9600 elements. 

The aforementioned optical calibration target provided a means to map the physical 
dimensions of the animal to the digital structure obtained from surface extraction. 
Specifically, the unit length in the CCD images was obtained by dividing the distance of the 
grids by the number of corresponding pixels. The coordinate system of the optical data was 
subsequently defined. Using the calibration procedure on the x-ray system described in [32], 
the unit length in the x-ray micro-CT data was determined. 

The x-ray micro-CT images were manually segmented into seven different tissue types: 
muscle, bones, lungs, heart, liver, stomach, and spleen using Amira software (Visage 
Imaging, San Diego, CA). Due to the way the animal was mounted on the rotation stage, large 
skin flaps were created, particularly between the forelimbs and the hindlimbs. These skin 
flaps were also identified in the segmented micro-CT data, and later excluded from optical 
sampling to avoid sampling and modeling problems. The segmented micro-CT image, which 
contains structural a priori information of the animal, was subsequently registered to the 
optically extracted surface using a lab-developed MATLAB (The MathWorks, Natick, MA) 
program by firstly scaling the data matrices to the same spatial dimensions, and secondly 
aligning external structural features of the animal visible in both imaging modalities. The 
segmented and registered micro-CT data served as the structural basis of optical 
reconstruction. The optical properties, i.e., the absorption coefficient ju a , the reduced scattering 
coefficient ju and the index of refraction n of different tissue types were adopted from the 
literature [33] and used in the forward problem solver (see Table 1). Note that the average 
values of the optical properties at the excitation and emission wavelengths were used. 



Table 1. Optical properties used in the Monte Carlo simulations: absorption coefficient 
(ji a ), reduced scattering coefficient (// '), and index of refraction (n) 





Muscle 


Bone 


Lung 


Heart 


Liver 


Spleen 


Kidney 


Stomach 


H a (mm -1 ) 


0.0349 


0.0242 


0.0672 


0.0275 


0.1623 


0.1623 


0.0311 


0.0069 


ju ' (mm -1 ) 


0.3709 


2.2929 


2.1040 


0.8875 


0.6371 


0.6371 


2.0661 


1.3560 


n 


1.37 


1.37 


1.37 


1.37 


1.37 


1.37 


1.37 


1.37 



Spatial registration of the x-ray data to the optical data enabled transformation of the x-ray 
coordinate system to the optical system, in which FDOT reconstruction was realized. Because 
the center of the field of view in the optical coordinate system lies on the rotation axis of the 
animal, the origin of the optical coordinate system can be related to the physical dimensions 
of the optical imaging platform. The origin of the source rays was defined as the spot at the 
center of the second mirror of the galvo-scanner; and the origin of the detector rays as the 
center of the CCD chip of the camera. These two coordinates were known by physically 
measuring the dimensions of the optical system (tolerance <1 mm). Knowing the origin of the 
source ray, the point where the source ray intersects the calibration target with known 
location, and the 3-D animal surface, the coordinates of the source that impinges the animal 
were computed (Fig. 1). Similarly, as soon as the detector positions were determined on the 
CCD images, the coordinates of the detector positions were also computed. 
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Fig. 1. Source localization using the calibration target. The rotation axis of the animal, i.e., the 
blue dashed-line in the 3-D view (left) and the blue dot in the 2-D view (right), is co-planar 
with the calibration target. The intersection of the laser ray with the calibration target (point a 
in the 2-D view) is known from the calibration image, which in turn determines the intersection 
point of the laser with the surface of the animal (point b). Note that the position of the 
calibration target and the origin of the laser ray are known. 

In free-space animal experiments, the air-medium interface is the animal skin in its natural 
state. As mentioned previously, due to the posture of the animal, large skin flaps typically 
were formed between the forelimbs and hindlimbs. These thin skin flaps, as well as other 
irregularities of the skin, caused significant problems for modeling the forward problem 
correctly, even using our Monte Carlo method. We implemented an anatomically guided 
sampling strategy to selectively exclude problematic sampling points from image 
reconstruction, illustrated in Fig. 2. These problematic sampling points were identified 
according to the following criteria: (a) a sampling point (either a source or a detector) was 
intersecting a skin flap; (b) any pixel within 0.5 mm (i.e., 11 pixels in our setup) from a 
sampling point in the emission or excitation CCD images was saturated; or, (c) the signal-to- 
noise ratio (SNR) at a sampling point was less than 2. The skin flaps were identified from the 
micro-CT images manually. The saturation of a pixel was characterized by a value of 4095 in 
the CCD image resulting from an all-one representation in the 12-bit analog-to-digital 
converter. The SNR at a sampling point was defined as the ratio of the mean to the standard 
deviation of the values within the sampling area. 

2.4. Image reconstruction 

2.4.1. Forward problem 

Image reconstruction in FDOT involves solving the forward and the inverse problems, both of 
which play important roles in overall quality of the reconstruction. As the forward problem 
solver, we used our Monte Carlo method previously reported in [15]. The Monte Carlo 
methods are a class of computational algorithms that obtain definitive results from repeated 
random sampling of physical or mathematical systems. The reliability of the Monte Carlo 
results is strongly dependent on the number of repetitions. From our previous experience of 
solving the forward problem in FDOT, the minimal acceptable number of photons simulated 
was on the order of 10 7 , which took -15 minutes to complete for each illumination source 
point using a single processor on a desktop personal computer (Dell Optiplex 745, Intel Core 
2 Duo 2.4 GHz processor, and 8 GB memory). In this study, we used a total of 240 
illumination points, each having 10 8 simulated photons, which would have required 600 hours 
to finish on our desktop personal computer. Another limitation of our Monte Carlo method is 
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the computer storage required for intermediate data that pass from the simulations to the 
construction of the Jacobian matrix. Our previously developed software package consisted of 
a Monte Carlo simulation module (implemented in C programming language) and a module 
for constructing the Jacobian matrix (implemented in MATLAB). Using 10 8 simulated 
photons in a typical geometry for the mouse, -1 gigabytes (GB) of data is produced for each 
illumination source. 




Fig. 2. Anatomically guided spatial sampling strategy based on structural a priori information 
from x-ray micro-CT. A representative pattern of spatial sampling (a grid of green dots) is 
superimposed on the surface (left) and the skeletons (middle) of the animal (based on micro- 
CT). On the right, the constituents of the segmented animal image are: (1) muscle, (2) bones, 
(3) lungs, (4) heart, (5) liver, and (6) skin flaps, where were used as the a priori information in 
the forward model; the spatial relations of the excitation laser (red arrow) and the emission 
fluorescence (green) with respect to the animal is illustrated, in which any emission or 
excitation position may coincide with the skin flaps and results in difficulty in reconstruction. 

Despite the above limitations, the Monte Carlo methods are the most accurate compared to 
alternatives, such as analytical approximation or numerical finite-element methods. To 
overcome the limitations of our Monte Carlo software, we unified the two separate computing 
modules into a single C language implementation which pass data between the two 
computational phases in main memory and thereby eliminate the need for intermediate disk 
files. The integrated software package was subsequently re- written to use the message passing 
interface (MPI) programming library in order to run on a large-scale parallel supercomputer. 
The newly implemented parallel Monte Carlo software was validated by confirming 
replication of the same results as the earlier implementation. 

Parallel computing was performed on an SGI Altix 4700 shared-memory system, which 
consists of 36 computer "blades." Each blade houses a pair of Itanium-2 Montvale 9130M 
dual-core processors and the four cores on a blade share 8 GB of local memory. A maximum 
of 132 cores can be allocated for a single run. The processors are connected by a high-speed 
"NUMAlink" interconnect, by which the local memory on each processor is accessible to all 
of the other processors. Each processor runs an enhanced version of the SuSE Linux operating 
system. The reduction of computing time scaled linearly with the number of allocated cores. 

2.4.2. Inverse problem 

The mathematical formulation modeling FDOT is very similar to that for perturbation DOT 
[4,9]. Therefore, the same method can be applied to compute the Jacobian matrices for both 
imaging modalities. In this study, we used our Monte Carlo algorithm previously developed 
for DOT to compute the Jacobian [34] : 
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where x is the quantum yield of the fluorophore; Ufi and U ex are the fluence rate measured at 
the fluorescence and excitation wavelengths, respectively; K is the complex wave vector; L p>m 
is the total path length of photon p in medium m; [F] r is the value of function F evaluated at 
location r; N s is the total number of simulated photons; and y is a constant scaling factor that 
accounts for the hardware- specific characteristics, such as the insertion loss of the filters, gain 
of the camera, and the optical power of the light source. 

The fluorescent image was reconstructed by solving the inverse problem using our 
previously developed regularized algorithm [21]. Image reconstruction was achieved by 
iteratively solving an unconstrained minimization problem 

X =argmin[#(X)] (2) 

where X is the fluorescent image, and H is a cost function regularizing the reconstruction. The 
cost function H(X) is defined as 

H(X) = \\AX-bf+a[ \xfdr + p[ \VXfdr (3) 
v 7 [] [5 Jq" 112 Jq" 112 

data consistency intensity penalty smoothness penalty 

where llvll 2 is the L2 norm of vector v, Q is the region of interest, and r is location. The three 
terms of the cost function penalize the residual error (i.e., data consistency), intensity, and 
smoothness of the solution, respectively, whose weights are independently adjusted by the 
user-specified parameters a and /?. This minimization problem was solved using the conjugate 
gradient method with the reconstruction space constrained within the torso of the animal. 

3. Results 

The matrix size for Monte Carlo simulation (the forward problem) was 186 x 135 x 260 with 
an isotropic voxel size of 0.18 mm. The resulting Jacobian matrix for the inverse problem was 
down-sampled by a factor of 4 in each dimension, producing a matrix size of 46 x 33 x 65 
and an isotropic voxel size of 0.72 mm. As a result, the Jacobian matrix occupied 7.6 GB of 
computer memory for 64-bit data precision. Using 128 cores of the supercomputer, the 
forward problem of the entire data set can be solved in -200 minutes. Overall, the 
computation performance scaled linearly with the number of computing cores. The inverse 
problem solver was implemented in MATLAB and took -15 minutes to finish on our desktop 
computer workstation (Dell Precision 690, dual Intel Xeon E5320 1.86 GHz 64-bit 
processors, and 40 GB memory). The inverse problem solver typically allocated up to 20 GB 
of memory. The regularization parameters, a and P in Eq. (3) were both empirically 
determined to be 0.1. 

The primary goal of this study was to demonstrate the capability of the proposed method 
to resolve two closely positioned small QD inclusions in the mouse thorax. The QD 
inclusions were visible in the micro-CT image, which measured 0.88 mm in diameter, and in 
length, 1.6 mm (for the rostral inclusion) and 1.8 mm (for the caudal inclusion). Figure 3 
depicts the axial slices of the FDOT reconstruction (normalized to unity) compared to the co- 
registered micro-CT slices. The slice thicknesses of the slices are 0.72 mm in FDOT and 
0.088 mm in micro-CT. The inter-slice distance is 0.72 mm in both imaging modalities. 
Although an imaging artifact exists in reconstruction, the two closely positioned microliter- 
sized QD inclusions were successfully resolved in the FDOT reconstruction. 
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Fig. 3. FDOT reconstruction (the first and the third rows from the top) compared to the 
corresponding x-ray micro-CT slices (the second and the forth rows). The FDOT data was 
normalized to unity, in which the small negative values were due to reconstruction artifacts. 
The actual inclusions were visible in the micro-CT image (rendered in red). The inter-slice 
distance is 0.72 mm. 

The localization error, defined as the distance from the centroid of the reconstructed inclusion 
to its actual location, measured 0.6 and 2.2 mm, respectively. We theorize the localization 
error was mainly due to the limitation of the forward model, as a result of the discrepancy 
between the assumed and actual optical properties of tissue. In principle, it is possible to 
determine the actual optical properties of tissues using DOT, which in turn can be used as 
functional a priori information for FDOT reconstruction, e.g., using the techniques described 
in [35-37]. However, to the best of our knowledge, a high-resolution DOT reconstruction 
technique has not been previously reported for animal experiments. Instead, we assigned 
literature values to the high-resolution anatomical structures obtained from micro-CT. In 
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addition, the organs of the animal were assumed to be homogeneous and clearly separated in 
the forward model, whereas in reality, not only are the organs highly heterogeneous, 
particularly in the lungs, but also the tissue boundaries were often ambiguous in image 
segmentation. 

We define the resolution in a similar way to the Rayleigh criterion: the maximum of one 
peak coincides with the half-maximum of the other peak. If these two peaks are Gaussian- 
shaped with equal amplitude, this definition translates to a dip of 15% between the two 
maxima. Based on this definition, we used a threshold of 85% to evaluate the shape metrics 
by measuring the diameter and the length of the reconstructed inclusions, shown in Table 2. 
Note that the dimensions of the inclusions are comparable to that of the voxel size, and that 
each inclusion in perfect reconstruction should occupy only 2 voxels. 



Table 2. Shape metrics of reconstructed QD inclusions 



Inclusion 


FDOT 


Micro-CT 


Diameter (mm) 


Height (mm) 


Diameter (mm) 


Height (mm) 


#1 (rostral) 


0.72 (1 pixel) 


1.44 (2 pixels) 


0.88 


1.6 


#2 (caudal) 


1.44 (2 pixels) 


1.44 (2 pixels) 


0.88 


1.8 



4. Discussion and Conclusions 

Although parallel Monte Carlo implementations have been previously reported by several 
groups, e.g., in [38,39], our parallel implementation is based on a different algorithm to 
construct the Jacobian matrix. As a result, our implementation is computationally more 
efficient when the number of detectors is large because it does not compute the photon density 
field of the detector. In addition by eliminating the intermediate data, applying our Monte 
Carlo method to preclinical FDOT with high-density optical sampling became practical. In 
conventional Monte Carlo implementations where photon density fields of both source and 
detector are computed, it is implicitly assumed that the source and the detector are 
interchangeable. This assumption is questionable for free-space FDOT because the source is 
typically a laser, which is highly directional and extends several millimeters into the medium 
(evident in Monte Carlo simulations, data not shown); and the detector is a surface area on the 
medium, which is much less directional (with a large acceptance angle). This assumption is 
particularly problematic for small animal imaging (e.g., in the mouse) where the 
source/detector distance is only an order of magnitude larger than the transport mean free 
pathlength (10-20 mm vs. 1-2 mm). 

Monte Carlo methods are inherently parallel because of the highly independent nature of 
simulation processes. In our implementation, the entire simulation consists of multiple 
acquisition angles, which was specified by the user. From the perspective of simulation 
software, each angle was a self-contained simulation job, which consisted of a number of 
sources and detectors. The simulation software divided the job according to the number of 
sources. For each source, the computer evenly distributed the workload (number of photons) 
onto the allocated computer nodes. Because the number of photons was orders of magnitude 
larger than the number of computer nodes (e.g., 100 million vs. -100, as in our case), the 
computing workload was almost perfectly balanced, i.e., all the nodes finished their tasks 
virtually simultaneously. Upon completion of the current source, the results were saved to the 
corresponding output file; and thereafter simulation for the next source started. 

The optical properties of biological tissue reported in the literature are highly inconsistent, 
due to the differences in experimental conditions and measurement methods. They are 
particularly scarce for near infrared. One should also recognize that the differences among 
species, individuals, and physiological conditions can be significant. Although the optical 
properties that we referred to in this study were for bioluminescence, we found that their 
model and parameters produced reasonable values in near infrared compared to those 
available in the literature. We decided to use the referred model and parameters from the same 
source for consistency. 
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In this study, we used the same averaged values of optical properties for both the 
excitation and emission wavelengths using our Monte Carlo software. Note that in 
conventional Monte Carlo methods, the Jacobian matrix was constructed by taking the 
product of the fields of photon density waves irradiated from the sources and the detectors at 
the excitation and emission wavelengths, respectively. In our preliminary comparison, we 
found that our method produced better image quality (data not shown). We hypothesized that 
it was because our Monte Carlo algorithm more closely represented the physical process of 
photon migration from the sources to the detectors within the medium. Further study is 
underway to better investigate our findings. 

Resolution is the ability of an imaging system to resolve two distinct objects, which is 
highly subjective to the viewer. In this study, resolution was defined as the center of a 
Gaussian peak coincided with the half-maximum width of another identical Gaussian peak, 
which translated to a threshold hold values of 85%. Note that different criteria and peak 
profiles can result in more or less different definition of resolution, e.g., 80% threshold if 
using the Rayleigh criterion and 95% using Dawes' limit for Airy diffraction patterns. 

One limitation of this FDOT study was the reconstruction artifact shown in Fig. 3. This 
artifact was likely due to the inaccurately known in vivo optical properties, especially in the 
lung area where the background signal was relatively high. Nonetheless, in preclinical 
applications, the disease model is typically known, which assists the interpretation of imaging 
results and avoids misjudgment. Note that under comparatively more extreme experimental 
conditions, such as the case in this study, imaging artifact tend to manifest more severely. It is 
certainly a direction for us to continuously work on to improve reconstruction quality in 
future. 

Another limitation was the lack of quantitative results. Although in this study, it was 
feasible to calibrate the scaling factor y in Eq. (1) using another QD solution with known 
concentration, we noticed that this calibration technique could not be readily generalized to in 
vivo animal experiments, where this calibration technique would become impractical because 
creating a fluorescent inclusion in the thorax of the animal is difficult, and is also problematic 
for preclinical studies. Another concern for this calibration technique was the use of the 
regularization method in reconstruction due to the ill-posed inverse problem. In addition, it is 
known that regularization in reconstruction affects the quantitative result in diffuse optical 
imaging [40]. The values of the regularization parameters (0.1 for both) indicated that the 
reconstruction was non-negligibly influenced by the second and third terms of the cost 
function defined in Eq. (3). With improved experimental control, signal acquisition, and 
imaging apparatus, we anticipate a decreased noise level in the optical signal and, as a result, 
significantly less regularization will be required in the reconstruction. In addition, more 
accurate optical properties will likely reduce the systematic error in the forward model. Under 
these circumstances, it may be feasible to use this calibration technique to produce reliable 
quantitative FDOT reconstruction for in vivo studies. 

In conclusion, we demonstrated a new FDOT method consisting of our newly developed 
parallel Monte Carlo software and anatomically guided sampling strategy for preclinical 
imaging. The Monte Carlo software fully utilized the power of parallel supercomputing and 
thus is better suited for handing large data sets in preclinical imaging. The anatomically 
guided sampling strategy significantly improved fidelity of data acquisition and subsequently 
the quality of reconstruction. Using the proposed method, the new FDOT system achieved 
significantly better resolution and sensitivity compared to our previous system, which 
consisted of essentially the same hardware. This method can be applied to any FDOT system 
adopting similar imaging principles to improve resolution and sensitivity. 
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